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513 ■ Abstract 



We present a generalized formulation of the Optimized Effective Potential (OEP) 
approach to the Self Interaction Correction (SIC) problem in Time Dependent (TD) 
Density Functional Theory (DFT). The formulation relies on the introduction of a 
^ ' double set of single electron orbitals. It allows the derivation of a generalized Slater 
^^ . approximation to the full OEP formulation, which extends the domain of validity 
^vq \ of the standard Slater approximation. We discuss both formal aspects and practical 
Tij" ' applications of the new formalism and give illustrations in cluster and molecules. The 
f--. ■ new formalism provides a valuable ansatz to more elaborate (and computationally 
^^ \ very demanding) full TD OEP and full TD SIC calculations especially in the linear 
I ' domain. 
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1 Introduction 



Density Functional Tlieory (DFT) lias evolved to be one of tlie most powerful 
theoretical frameworks for the description of complex chemical and physical 
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systems. Enormous progress has been made since the seminal works of the 
sixties by Kohn et ah [Il[2]. DFT is now routinely used, especially in systems 
with a large number of electrons [3HS]. Nevertheless, there remain still several 
open questions in detail which are in the focus of actual research j^]. The 
extension to Time-Dependent DFT (TDDFT) has been formally established 
more recently [3-[9]. It is still a developing field, at the side of both formal 
and practical aspects [10]. Already at the present stage, TDDFT has become 
one of the few well founded theories for describing the dynamics of complex 
systems. This is especially true for non equilibrium situations such as clusters 
and molecules under the influence of intense laser fields [TT] . 

DFT simplifies the involved problem of many-electron correlations in terms 
of an effective (Hartree-like) one-body description. This is achieved by intro- 
ducing exchange and correlation effects in an energy functional expressed in 
terms of the local density of the electrons. The simplest strategy along that 
line is provided by the Local Density Approximation (LDA) which has proven 
in many calculations to provide a simple and reliable description of structure 
and low-amplitude excitations (optical response, one-photon processes) [5]. 
The analogue in the time- dependent case is the Adiabatic Local Density Ap- 
proximation (ALDA) which has also been used with great success in dynamical 
processes involving huge energy deposits and/or large ionization of irradiated 
clusters and molecules [TT] . 

However, the LDA is plagued by a self-interaction error because the direct 
Coulomb term and the LDA exchange-correlation potential involve the total 
density including the particle on which the mean-field acts in the Kohn-Sham 
equations. While in a full Hartree-Fock treatment, the exchange term exactly 
cancels the self interaction of the direct term, the approximate treatment of 
the exchange term in LDA destroys this cancellation. As a consequence, a 
self-interaction remains and one of the defects is that LDA produces a wrong 
Coulomb asymptotics [11[I2]- The self-interaction thus spoils single-particle 
properties in particular the Ionization Potential (IP) in finite systems or the 
band gap in solids [T3|[T^. It is also well known that LDA fails in describing 
the polarizability in chain molecules [T5l[T6]. In the dynamical case, the self- 
interaction error also spoils the description of ionization dynamics, especially 
close to threshold where IP effects dominate. 

There exist ways to correct the self-interaction error while trying to keep the 
simplicity of the method. An early attempt along that line was proposed by 
Fermi and Amaldi [T7] . The standard way to introduce a Self Interaction Cor- 
rection (SIC) is based on the more recent proposal of Perdew [T^fTS] . Such 
SIC has been explored since then at various levels of refinement for struc- 
ture calculations in atomic, molecular, cluster and solid state physics, see 
e.g. [T9] - [22] . The SIC scheme, however, leads to an orbital dependent mean- 
field which causes several formal and technical difficulties. A way out is pro- 



vided by using optimized effective potentials (OEP) techniques as introduced 
in [23l|23], see [25] for a recent review. However, applying OEP to SIC leads 
to a very involved formalism usually treated with further approximations, as 
e.g. the Krieger-Li-Iafrate (KLI) approach [2^127] . But these approximations 
can severely perturb some crucial physical features of SIC, particularly the 
trend to produce localized single-particle orbitals [25]. It is thus a key issue 
to refine such approximate schemes to SIC in order to preserve, as much as 
possible, original SIC properties, at a lower cost than full OEP. It should be 
noted that there is nevertheless a further advantage of OEP. It optimizes one 
local mean-field Hamiltonian for the system. This allows to evaluate unam- 
biguously unoccupied states of the system which, in turn, can have important 
applications in dynamical processes. 

Time-dependent situations call for a time-dependent SIC (TDSIC). Applica- 
tions of TDSIC have, up to now, mostly been performed in approximate man- 
ners, e.g., the linearized treatment of [2S], averaged-density SIC [22] based on 
a generalization of the Amaldi picture [17], or the various versions of time- 
dependent OEP-KLI [301432] . Only recently, a manageable and exact propa- 
gation scheme for TDSIC has been formulated [33l[3l] which is applicable in 
all dynamical ranges. The key to success is to employ two complementing sets 
of occupied single particle wave functions, one for the mean-field propagation 
and the other one establishing the necessary localization of the wave functions. 
The double-set technique has also proven to be extremely useful to formulate 
efficient approximations to OEP. This was demonstrated for the stationary 
case in [SHIES] and later on extended to the time domain [37]. The aim of the 
present paper is to present and discuss in more detail the local approximations 
to time dependent OEP (TDOEP) based on the double-set technique. 

The paper is organized as follows. In section [2l we briefly review the double- 
set technique for SIC and TDSIC. In section [3], we introduce the (TD)OEP 
equations in the light of the double-set representation and develop from that 
what we call the generalized Slater approximation to OEP. In section HI we 
present results for a variety of test cases, static as well as dynamic ones, and 
compare with results from the lower approximations LDA and ADSIC, and 
from full TDSIC as a benchmark. Conclusions are summarized in section [51 



2 Brief review of the double-set technique for SIC 



In this section, we give a brief outline of SIC and TDSIC in the double-set 
formulation as introduced in [331 |3l]. We start from the static case which 
helps to motivate the double-set technique and proceed to the dynamical case 
where the use of two sets of wave functions is compulsory. The brief review 
of (TD)SIC should serve as a starting point for the derivation of improved 



approximations to OEP (both in stationary and time-dependent cases) 

2.1 Stationary case 

The starting point is the SIC energy functional for electrons : 



whereby all sums run over occupied states only. Note that we omit the space- 
time dependencies, i.e. ipa = '4'a{^i'^)i when it is not misleading, p stands for 
the total electronic density p = J2aPa = Z^a I'^oP- The first term in E'sic 
is the non-interacting kinetic energy; E^xt [p] = I d^^pvcxt collects all external 
one-body fields where Vext stands for the interaction with the ionic background 
and any other local possibly time-dependent external field; and -EldaIp] is 
a standard LDA energy- density functional including the direct term of the 
electron-electron Coulomb interaction. The last term corresponds to the SIC. 
We mention in passing that the SIC, and with it all our following development, 
does also apply to more general functionals as, e.g., the Generalized Gradient 
Approximation (GGA) |38]. A basic assumption beyond this SIC functional 
is that the employed single-particle wave functions are ortho-normalized 

ii^al^Pp) = S^p . (2) 



The stationary SIC equations are obtained by variation of the SIC energy 
dl]) with the additional constraint on ortho-normalization (|5]) which is taken 
into account through the Lagrange multipliers Xa/s- This leads to the following 
mean-field equations p ^ [33 | IM] 

hsicli^a) = Yl ^apli^a) ) (3a) 

hsic = hLDA - Usic , (3b) 

p 

hhBA = - h UldaIp] , (3c) 

2m 

^lda[p]= ^^^^^ , (3d) 

f>SIC = Ef^LDA[|^ar]|^a)(V^a| , (3e) 

a 

combined with what we call the "symmetry condition" 

0=(V^/3|t/LDA[|^/3r]-f^LDA[|^a|']|V^a) ■ (3f) 



The first term in /isic is the standard LDA mean field Hamiltonian (l3c|) and 
the second term stems from the SIC in the energy ([T]). 

Thus far, Eqs. ([3]) comprise the complete stationary SIC method. The right- 
hand side of Eq. f pa|) is unconventional and inconvenient as it does not lead 
explicitly to single-particle energies. The latter may be defined a posteriori 
by diagonalizing the matrix of Lagrange multipliers X^is- We now introduce 
explicitly a second set of wave functions (Iv^i)} which indeed diagonalizes the 
SIC Hamiltonian, 



hsic\<^i) = £i\<-Pi) , (4) 

and which is related to the previous set of {|'?/'a)} by a unitary transform 
amongst occupied states only : 



l/ja =J2^i Mia , Yl ^ic^i^ = ^a/3 ■ (5) 

i i 

Both sets lead to the same total density p such that the LDA mean-field 
UldaIp] remains the same. The new set {Iv^j)} represents the energy diagonal 
states, while the old set (iV'a)} remains the decisive ingredient in the symmetry 
condition fl3fj) . The coefficients Uia of the unitary transformation ([5]) for given 
ipi are to be determined such that the symmetry condition fl3fj) . involving 
the ipa, is fulfilled. As the ipi orbitals satisfy eigenvalue equations, they are 
interpreted as single electron orbitals. The set ipa serves to minimize the SIC 
energy ([1]) and to calculate the SIC mean-field /isic- 

This completes the double-set representation of stationary SIC. The double-set 
technique is not compulsory for the stationary case, but enlightening. The two 
sets play different roles. The energy diagonal states can easily be delocalized 
and are likely to spread over the whole system, e.g., when considering the 
valence shell of metallic bonds. The SIC set {|^«)}, on the other hand, aims 
to minimize the SIC energy which is usually achieved by localization of the 
associated density l^ap to minimize the Coulomb energy 



2.2 Time- dependent case 



In contrast to the static case, the double-set technique is a necessary ingredient 
for developing a well-defined and manageable propagation scheme. 

To derive the TDSIC equations, we start from the SIC quantum action 



^sic = I ' dt(^Esic[{^Pa}]{t) - J2iMtMdt\Mt))) ■ (6) 

The situation with the action for time-dependent variation is, in fact, not so 
trivial for the derivation of DFT. There can arise problems with causality [ID] 
and boundary conditions [H] for which solutions are discussed in [HIHS] . We 
are dealing here with a local and instantaneous ALDA functional which allows 
to use the naive action ([6]). Moreover, concerning the theorems derived from 
symmetries of the action, it is shown in [12] that, as compensations occur, 
the stationarity of the naive action ([6]) leads to the correct final results. We 
thus perform variation of this action including once again the ortho-normality 
constraint with Lagrange multipliers A^/3, i.e. we require 



6(Asic- dtJ2{Mt)\Mt))\fs{t) 

It is to be noted that, to derive the time-dependent OEP formalism, one should 
use the action ([6]) in the limit to ~^ —00. This is necessary to recover in the 
stationary limit the stationary OEP formalism, as proved in ^3] . 



The steps of the variation are explained in detail elsewhere [33J[3lj. We sum- 
marize the resulting equations. They again employ the two sets of occupied 
single-particle wave functions which are connected by a unitary transforma- 
tion ([5]). The set {\(pi{t))}, which was the diagonal set in the static case, now 
turns out to be the "propagating set" obeying the time-dependent mean-field 
equation 



hsicit) - ihdt]\^i{t)) = , (7) 



where /isic is defined by ( l3ell . The coefficients Uia of the unitary transform 
([5]) for given ipi are to be determined such that the "localizing set" {{ipait))} 
satisfies the symmetry condition 



U^t) : Vt, = (^/3(t)|f/LDA[|^/3|'](t) - Ui^DA[\i'a\^]{t)\Mt)) (§) 

at any time. The solution scheme for these two coupled equations is obvious. 
The time- dependent Schrodinger equation ([7]) for the propagating set is solved 
for a short time step by standard techniques, e.g., a Taylor expansion of the 
formal solution \(pi{t')) = expj— ^// dr /isic(''") f|v^j(^))- At each time step, 
the set ipa is determined by resolving the symmetry condition (l3fjl [M] , which 
is an instantaneous equation. Then the ipa serve to construct the new mean- 
field hsic for the next time step. 



This TDSIC propagation scheme looks formally straightforward. However, 
it contains one especially numerical expensive ingredient: The iteration of 
the symmetry condition (l3fj) requires to invoke the time-consuming Coulomb 
solver in each iteration step. Any acceptable approximate solutions are thus 
welcome. The time-dependent generalized Slater approximation which will be 
discussed below is a step into this direction. 



3 SIC-OEP and the Generalized SIC-Slater approximation 

3. 1 Stationary formalism 

3.1.1 SIC-OEP in double-set representation 

The "Optimized Effective Potential" (OEP) formalism is the tool of choice to 
find the best local approximation to a non-local Hamiltonian. In the present 
case in which we plan to apply OEP to the SIC problem, we start from the 
total SIC energy ([1]) formulated in terms of the (localized) ipa orbitals and we 
complement this set by the diagonal orbitals (fi. The latter are required from 
the onset to satisfy a local eigenvalue equation 

[hLDA-Ui^^^\r)jipi{r)=eiipi{r) . (9) 

Locality is imposed by the fact that Uq^iq'^ (r) is a function of r only. The 
localizing set of tpa is obtained by the unitary transformation ([5]) whose coef- 
ficients are optimized to minimize the SIC energy ([1]). It remains to determine 
the space of occupied single-particle states in terms of the (pi. The condition 
(jH]) shifts the problem to a yet unknown optimizing local potential U^iq^ ■ 
This potential then becomes the variational degree of freedom instead of the 
Lfi. The potential U^iq^^ (r) is thus determined by minimization of the total 
SIC energy ([1]) 

^^^Jf#^ = . (10) 

Note that no additional ortho-normality constraint is needed in the variation 
because it is aheady guaranteed by solving Eq. (^. 

The variation is performed using the chain rule for functional derivatives 
5Esic/SUg°Q^ = {6Esic/S^i) i^^i / ^^sic^ ) where the first factor represents 
the usual SIC mean-field and the second factor the wave function response to 
varied local potential. The detailed derivation is given in [55] . We obtain as a 



final result an integral equation for the t/gi'c '^ 

E/dr'(f/Sr'^(r')-<(r'))G.(r,r')^:(r>.(r)=0 (11a) 

where 

Gi(T,Y) = }_^{l-5ij)^ — Z (lib) 

is the single-particle Green function in the mean-field (jH]). The driving quantity 
in the integral equation (11 lap is the SIC potential with respect to the (^j, 
namely 



1 ^^icKV^j] 1 ^ 

V9i(r) 5(^*(r) v5i(r) 

E<f^LDA[|V^.|'](r)V'a(r) = ^(r|f/sic|^.) , (He) 



where the ipa are deduced from the ^pi by the unitary transformation ([5]) 
whose coefficients are determined by the symmetry condition fl3fl) . Note that 
we considered in this variation the diagonal basis states (pi and the coefficients 
u of the unitary transformation to the ipa as independent, i.e. 6u/6(f*{r) = 
(see also section 4.4 of Ref. [M] where we have shown that the (fi and the 
Uia should be considered as independent in the variation of the energy or the 
action) . 

Eq. (II lap defines U^ic'^ (r) in a rather involved way. Its solution can be dis- 
entangled as 



f^slr'^ = ^s + V^K + V^c , (12a) 

^s = E^^^ ' (12b) 

i ' 

I |2 
^K = E —^i'^ilUilc'^^^ - Vi\ipi) , (12c) 

P.(r) = ^/dr'(4r'^(r') -<(r'))^*(rOG,(r,rO . (12e) 

From a practical point of view, this form is not simpler to use than the original 
form (II lap because the Vk and Vq terms depend on the solution U^jq^ . How- 
ever the separated representation serves as a starting point to develop further 
approximations. 



Some straightforward manipulations with the unitary transformation ([5]) allow 
to rewrite these quantities in terms of the localized wave functions ipa as 



a P 

^K = -E (E IV'^r<«^/5)(^/3|f^s!r'^ - ULDA[\i^ama) , (13a) 

P^(r) = ^E^- /dr'felr'VO - f/LDA[|^„P](r'))c(r')G.(r,r') . 

(13b) 

This expresses SIC-OEP in terms of the double-set representation what we 
call the "Generalized SIC-OEP" formalism. Thus far it is fully equivalent to 
the original SIC-OEP equations f lTT|) and as involved to solve. But the double- 
set form flT3|) employs the more localized states ^/'q which produces a more 
forgiving hierarchy of importance for the different terms. 

It is to be noted that a very similar development is found in [M] , but without 
addressing the feature of spatial localization of the V'a when introducing the 
KLI approximation. 



3.1.2 Generalized SIC-Slater approximation 

In this section, we show that the spatial localization of the ipa permits to 
justify a powerful approximation. We define 

12 
/3 P 



i^r^ = (u^^m'] - E ^^f/LDAiiv^.r])^^ . (14) 



In the expected case that the ifja remain spatially localized, we can assume 
that at each space point r one ifja dominates. This means that 

Fp) ^ . (15) 

We now take up the "Generalized SIC-OEP" equations ( 1T3l) and reshuffle them 
to display the FJ^^'^ explicitly. The Vg term is dominating compared to the 
Vk and Vq terms, although those latter terms may be not a priori negligible. 
Thus we assume approximately 

^sr'^^v^s = E^f^LDA[i^«r] • (16) 



Inserting Eq. ( !T6|) into Eqs. fll3ap and (113bp . we obtain for the two remaining 
pieces 

P^(r) = -^E^-/dr'i^P"*(r')G.(r,r') . (17) 

Using the feature f lTSj) which follows from the localization of the t\}a^ we obtain 

V^K^O , p, ^O^yc~0 . (18) 

This justifies a posteriori the assumption ( TT6|) and so allows to neglect the 
more involved contributions Vk and Vq [35]. Note that the double-set technique 
leaves full freedom for the diagonal orbitals (y9j, whose degree of localization 
can strongly vary according to the studied system (the (y9j are, e.g., strongly 
delocalized in a metal and more localized in covalent binding). 

After all, the generalized SIC Slater (GS) approximation to SIC-OEP [35] can 
be summarized in the three coupled equations 



^LDA -f/Gsj|v5i)=ei|v5i) , (19a) 

t/GS=E^^^LDA[|^aP] , (19b) 
a r 

= (^.|?7LDA[|^a|V^LDA[|^/3|']|^/3) • (I9c) 

Eq. (I19ap determines the (pi for given Uqq,- Eq- (I19cp determines the localized 

states ipa by finding an appropriate unitary transformation (j5]). These if^a are 

employed in Eq. fll9bp to determine the local mean- field Ugs- The coupled 
equations can be solved by iteration [35] . 



3.1.3 Comment on the traditional SIC-Slater and SIC-KLI approximation 



One can show that the traditional SIC-OEP formalism [261127] is obtained 
by the same reasoning as previously, i.e. imposing that the diagonal orbitals 
i-Pi satisfy a local Schrodinger-like equation (I9l). But having no second set of 
wave functions at hand, it stops at the stage (fT2l) . The traditional SIC-OEP 
equations [26l|27] are then obtained replacing (lllcp by 

Vi = Ui^BA[\^i?] (20) 



10 



in the equations (11 lap and (112bp - fll2ep . The traditional SIC-Slater approxima- 
tion is here formulated as 



I |2 
t/s = E^^^LDA[|¥^^r] (21) 

i r 

We discuss its validity in terms of the quantity 



fP = (u^daHv.]'] - E ^f/LDA[|^,f ])^. . (22) 

One can shows that the Slater approximation is justified only if 



fP^O , (23) 

i.e. only if the (fi orbitals remain spatially localized or very delocalized (close to 
a Fermi gas). Similar reasoning applies to the traditional KLI approximation. 
But as extensively discussed previously, the (pi have no particular reason to 
remain localized or to delocalized in the general case. There are many favorable 
situations as, e.g., a tendency to localized orbitals in organic molecules but 
also strong delocalization in metallic systems. Thus there is quite a choice of 
systems where the Slater approximation is found to be applicable. However, 
there is also a great number of systems where the Slater or KLI approximations 
fail. This has been numerically shown in |13] for the example of hydrogen 
chains. 

The GS approximation contains with the double-set technique an extra local- 
ization step which significantly enhances the range of validity of the Slater ap- 
proximation. This was shown in terms of several numerical examples in [35|[36] 
where, e.g., the demanding hydrogen chain was found to be reasonably well 
described within the GS approximation. A key issue for justifying a Slater- 
type approximation is that the single electron LDA term L^ldaIV'o] is close 
enough to the density weighted average (po/p) thereof. The approximation 
obviously works when the ipa are sufficiently localized, as, for then, around 
the given point where ipa is localized, one has p c:^ p^. These strongly local- 
ized orbitals correspond to a hydrogen or rare gas bond. There is the other 
extreme of metallic behavior in which all ipa extend over the whole system 
and whose densities resemble each other. This also provides an a priori well 
working Slater approximation. In between these two extremes range numerous 
conceivable cases. In particular, they can involve covalent binding which are 
thus not so well approximated in a simple minded Slater picture. We shall see 
below that the double set formulation allows to address also such intermediate 
bindings and performs well in these cases as well. 
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3.2 Time- dependent formalism 

3.2.1 Time- dependent Generalized SIC- OEP 

We now develop the time-dependent SIC-OEP and "Generalized SIC-Slater" 
formalisms [37]. Starting point is again the SIC quantum action (E]). We impose 
that the orbitals ipi satisfy a time-dependent Schrodinger-hke equation with 
local mean-field potential f/gj^^ (r,t), i.e. 

hi^BA{r,t)-Ui°^'''\r,t)-ihdty.{r,t) = . (24) 

The optimal f/gj^^ (r,t) is to be determined by variation of the action 

^^^^^ (25) 



6Ui-\r,t) 

while the single-particle wave functions ipi become potential-dependent quan- 
tities. No additional ortho-normalization constraint is needed in this variation 
because it is already provided by the solution of Eq. (IMl) . The localized set V'a 
is again deduced from the ipi by the unitary transformation ([5]) and the trans- 
formation coefficients are to be determined by variation of the action. This 
yields once more the symmetry condition fl3fj) to be fulfilled at each instant. 
It is to be noted that the emerging double set of ipi with ip^ is not exactly 
the same as the solution of TDSIC Nonetheless we use the same notations for 
sake of simplicity. 

Similarly as in section 13.1.11 the variation ( l25l) is again evaluated with the 
chain rule for functional derivatives. After a series of formal manipulations, 
one obtains an integral equation for the optimal local mean-field potential 

TT(local) 
'^SIC 

= E£jt'Jdr'[uiS^'\r\t')-v*ir',t')y<,{r,t;r\t')^:^^^ 

+C.C. , (26a) 

K,ir,t;r',t') = -t ^ ^*{v,t)^j{r',t')e{t - 1') , (26b) 

^.(r,t) = ^ .. , !, .. f dt'Esicit') + —1— (r|/iLDA(t)|(p.(t)) 

ipi{r, t) dv5*(r, t) J-oo ipi{v, t) 

E<W^LDA[|^a|'](r,t)^<.(r,t) . (26c) 



V5i(r,i) 



7-(local) 



As in the static case, we can decompose U^iq in terms of separate contribu- 
tions 
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f^slr'^ = ^s + ^e{VK + Vc} - 53m{V^TDi + 1^td2} (27a) 

where Vs, Vk, Vq are expressed exactly as in (112bp - (112dp but where now the 
time dependence induces possible complex components we shall analyze fur- 
ther below. However Eq. (I12ep defining the Pi for Eq. (112dp is to be replaced 

by 



pU,t) = -^-°^ ^ ^ 7 ^ (27b) 

The potentials Vs, l^, Vc, which also appear in the stationary case, are now 
complemented by two dynamical contributions 



^TDi = -E^^r dt\ip.{t')\v,{t')\^.{t')) , (27c) 

V^^, = -Y.{\v^?^-^ + h.Vp^ , (27d) 

where Jj = ^(v9*Vv9j — ipiVip*) is the current density. Note that the potential 
Vtdi contains a time integral, thus memory effects, while the potential Vtd2 
involves the time derivative of the pi. 

The standard way to derive the time-dependent Slater and time-dependent 
KLI approximations starts from the above separation in terms of the propa- 
gating basis ipi. More robust approximations will be obtained from a formu- 
lation in terms of the localizing set ipa- The separation can be remapped to 



V^s = E^^f^LDA[|V^«P] , (28a) 

a P 

^K = -E (E \V^\'nlu,,){MUi^^'^ - U^DAil^Pama) , (28b) 

P a,l3 ^ i ^ 



P.(M) = -i-EMt)/^dt'/dr' 



</'r(r,t) 
f/Sr'^(r',t')-f/LDA[|V^.P](r',t'))c(r',t')^.(r,t;r',t') 

(28c) 



-aPll^a) , 



{^i\Vi\^i) = E"i/3^*L(^/9|f^LDA[|^o 

= E"i""i'/3(^a|f^LDA[|^aP]|V'/3) = {^ilVilipi)* ■ (28d) 
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The last equation (128dp implies Qm{{ipi\vi\ipi)} = which, in turn, yields 
Sj?7i{Vtdi} = 0. This thus allows to remove this term in the decomposition 
(I27iD. 

These equations for the optimal local mean-field together with the time- 
dependent mean-field equation (!24l) and with the symmetry condition (l3fl) 
constitute TDSIC-OEP in double set representation. Its solution is by no 
means simpler than the solution of fully fledged TDSIC But the equations 
with explicit separation of the optimal local mean-fleld provide a good starting 
point for approximations. 



3.2.2 Time- dependent generalized SIC- Slater approximation 

The reasoning to derive a time-dependent generalized Slater approximation 
proceeds very similar to the static case (section [3.1.2p . We introduce the func- 
tion F^'^^\v,t) deflned in flT^ and assume the generalized Slater approxima- 
tion 



^Sr'^^V^s = E^f^LDA[|^„P] . (29) 

a P 

Inserting into Eqs. fl28b|) and fl28cl) yields 



Vk = --Y.(J:\^^\"<^^^ /drFf S)(r,t)^;(r,t) ^ 



p,(r,t) = -^^S2^Ut) f dt'/dr'Ffs)*(r',t')^,(r,t;r',t')^0, 
^ \/td2 ~ , Vc^Q . 

These results are consistent with the approximation ( l29l) . Thus the potential 
( l29l) is probably a good approximation of the time-dependent SIC-OEP po- 
tential for a broad class of problems, as for example when the localization 
by the symmetry condition ( l3fl) works well. The time- dependent generalized 
Slater potential ( 129|) has the same form as the stationary GS potential. It 
does not contain memory effects anymore, which is another consequence of 
the localization of the il)a. 

The emerging scheme is called time-dependent generalized SIC Slater approx- 
imation (TDGS). It can be summarized by the coupled equations 
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{hLDA-UGs]\^i)=i^dt\ipi) , (30a) 

lib P 

f/GS=E^^^LDA[|^.P] , (30b) 

a P 

Vt : = (V'a|f/LDA[|^ar]-f^LDA[|^/3|']|^/3) , (30c) 



where the ipa are obtained from the ipi by the unitary transformation (j5]) which 
makes the ipa to satisfy the symmetry condition (I30cp . 



For the same reasons as discussed in the static case, TDGS should represent 
an improvement to conventional time-dependent SIC-Slater and SIC-KLI ap- 
proximations to the extent that it opens a larger class of problems for which 
the approximation is applicable. This will be demonstrated on practical test 
cases in section HI 



3.2.3 Conservation law I: Energy conservation 



Within TDGS, the (pi orbitals propagate under the influence of the potential 
fl30bl) according to Eq. O30a|] . The total energy is computed with Esic as given 
in Eq. ([1]). We remind that variation of ii^sic defines the SIC mean-field Usic 
as defined in (l3eD : 



^ -Esic[{^j](t) = (r|/isic(t')l^.(0) ^it-t') 



5^t{v,t' 

hsic{t) = kDA{t)-Usic{t) ■ (31) 



Energy conservation is an issue for time-independent external fields, i.e. for 
dt^ext = 0. The time evolution of the energy thus becomes 



15 



5 



= Y. J ^i'{dt^iit')\hsic{t')^r{t')) 6{t - t') + c.c. 

i 
= \Y. {hGsit)^iit)\hldt)^i{t)) - (hslc{t)^^{t)\hGsit)^^{t)) 

+|^/drE^f^LDA[|^,P]E^^A^:} 






We finally obtain the time variation of the SIC energy ([T]) within a GS prop- 
agation 



a*i?sic[{^J] = --Sm|^ /drFf s)(r,t)AC(r,t)| (32) 

where we employed the deviation function F^^'^ from Eq. (1140 . The relation 
(1321) shows that the energy is not strictly conserved. The quality of energy 
conservation depends on the quality of the generalized Slater approximation 
because the deviation is driven by the same function F^^^^ which enters the 
decision on negligible terms at the end of section I3.2.2[ 

We compare this result with the traditional time-dependent SIC-Slater prop- 
agation. We can show similarly that the time variation of the associated total 
energy £'sic[{v'i}]; now expressed in terms of the diagonal orbitals ipi reads 



^tEslc[{v^}] = -^Sm|^|drFf (r,t)Av,*(r,t)| (33) 

with Fl given by Eq. ( 1221) . Energy conservation holds for the traditional 



(S) 

SIC-Slater scheme if -Fj ~ 0, thus if the physical system as a whole re- 
mains homogeneous or very localized. The extra localization for the ip^ in the 
double-set technique makes it very likely that energy conservation is improved 
for TDGS. It is known for traditional TDKLI and TD-Slater that energy con- 
servation lasts for a certain time interval after which energy explodes |15]. We 
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expect that TDGS has a similar behavior but with a much extended time in- 
terval of practical energy conservation, which allows to use it in a wider range 
of physical situations. 



3.2.4 Conservation law II: Zero Force Theorem 

The Zero Force Theorem (ZFT) states that a time variation of the total elec- 
tron momentum can be caused only by an "external" potential [ISHU], i.e. 



dt{P) = - J drpVv,,, . 

It stems from the fact that the electron-electron interaction is translational 
invariant and can not produce a "net" force on the system which, in turn, 
leads to the ZFT in the form [17] 

Vp : Jdrp{r,t)VUrn{[p]{r,t) = (34) 

where t/mf is the local mean-field potential of the considered method. The 
ZFT holds for the LDA mean-field and the ADSIC one. We now are going to 
check the ZFT for the TDGS mean- field t/cs- 

The time evolution of the total momentum is given by 

i i 

+ (r|i)cxt + f^LDA - t/Gs|v^j)Vv?* (r) > 
= I dri;cxtVp + / drt/LDA[p] Vp 

-/dr^^f/LDA[|V^ar]Vp (35) 

The second term disappears as can be shown by a partial integration and 
exploiting the ZFT for ?7lda- We add = Ea I drUi^DA[M'^]^\ipa\'^ to the 
third term and reshuffle Eq. (1351) to 

9iE(<^dPl<^^) = -/drpV^;ext-23f?e|5:y'drFfs)(r,t)VC(r,t)} (36) 

where we employ -F^*^^) from Eq. (fT4|) . Again, we see that the deviation func- 
tion F^*^^) drives also the term that violates the ZFT. The ZFT is well fulfilled 
if FJ^*^^^ is small, i.e. if TDGS is valid. In reverse, violation of ZFT and energy 
conservation is a valuable indicator for the breakdown of TDGS. 
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With similar steps we can evaluate the ZFT for the traditional time-dependent 
Slater approximation and obtain 



ai5:(<^,|p|<^,) = -/drpV^;ext-23fJe|^/drFf)(r,t)Vv^*(r,t)| (37) 

i i 

where we use the FJ; from Eq. (1221) . The ZFT is thus verified within a tra- 
ditional SIC-Slater propagation only if Ff ^ 0, thus if the physical system 
remains homogeneous or very localized. We have argued above that the range 
of -Fj f^ is much smaller than the range of F^'^^'^ ^ 0. This means that 
TDGS should maintain the ZFT for a longer time span than traditional time- 
dependent SlC-Slater. 



3.3 Alternative localization criteria 



One major effect of the symmetry condition (I30cl) is that it produces states 
i/jce which are better localized than the originally given ipi. This was the par- 
ticular feature which we employed to motivate TDGS. On the other hand, the 
symmetry condition is very expensive to solve in practical calculations. Thus 
it is worth trying to achieve better localization by less demanding equations. 
There exist many localization criteria [121120] • After a series of numerical tests 
with many of these localization criteria, we have found as a best compromise 
for a localization criterion the spatial variances of the one-body orbitals : 

A^sp = E [(V^«|r>„) - Mr\^a?] , (38) 



where the index "sp" stands for the summed single-particle variances. Mini- 
mization of this variance yields the localization equations 

= (^a|r^- rsl^/j) , r^ = {ipa\r\ipa) (39) 

which then replaces the symmetry condition (I30c|) in the TDGS equation. It 
serves to determine the coefficients of the unitary transformation i^. It is 
again a non-linear equation which has to be solved iteratively. But the expec- 
tation value r^ can be computed much faster than the Coulomb field. Thus 
TDGS with the localization condition ( 139|) is computationally less demanding. 
We have to see how it performs in practice. 



4 Numerical results 



4.1 Brief reminder of the various studied formalisms 



In the following, we will compare the results obtained with (TD)GS and other 
approaches to those obtained with full (TD)SIC as a benchmark. The corre- 
sponding mean-field Hamiltonians are summarized in tabled! all being used in 
one-body Schrodinger-like equations of the form h\(pi) = ihdt\(fi). Note that for 



expression of /i in /i (/?j) = ei\ipi) 


method 


acronym 


hhDAip] 


LDA 


LDA 


hhDA [p] - Ulda -tt 




Average Density SIC 


ADSIC 


/iLDAH->;^^LDA \V,\' 


Standard SIC Slater 


Slat 


= (^If/LDAil^P] - f/LDA[|^aP]|^a) 


Generalized SIC Slater 


GS(sym) 
GS(var) 


a 

= (^/3|^LDa[|#|'] - f/LDA[|^a|']|^a) 


full SIC (benchmark) 


(TD)SIC 



Table 1 

The hierarchy of mean-field Hamiltonians, from simple-most LDA (top line) to full 
TDSIC (bottom line). The right column shows the acronyms used in the figures and 
discussion. 



GS the symmetry condition ((Sfl) should be added for the two last schemes, to 
define the localized states ifja required in the corresponding Hamiltonians. As 
an alternative, we consider the localization criterion f l39|) derived from mini- 
mization of the spatial variance ( l38l) . see section l331 We abbreviate the scheme 
using the symmetry condition as "GS(sym)" and the alternative scheme using 
Eq. (EHD as "GS(var)". 

The static and dynamical calculations are performed on 3D coordinate-space 
grid using standard techniques, for details see e.g. [UlEI]. The calculations 
are restricted to valence electrons. They are the 3s electrons in Na, the 2s and 
2p electrons in C, and naturally the Is electron in H. The coupling of the ionic 
cores to the valence electrons is described by pseudopotentials. For the C and H 
atoms, we use Goedecker-type pseudopotentials [52] and for Na atoms the soft 
local pseudopotentials of [53] . The LDA part employs the exchange-correlation 



19 



energy functional from [5l]. ADSIC is performed as explained in [29]. The 
static solution is done by accelerated gradient iteration [55|[56]. Time stepping 
is done by fourth order Taylor expansion of the exponential propagation oper- 
ator [57] . The Poisson equation is solved by a fast Fourier technique combined 
with separate treatment of the long-range terms [58]. Polarizabilities are com- 
puted from two static calculations where one is performed under the influence 
of a small static external dipole field. 



4-2 Potential energy surfaces 



The C2 molecule is found to be a critical test case. The electronic structure 
changes substantially from the spin-saturated, covalently bound dimer ground 
state to the highly spin polarized asymptotic atomic states. It is demanding 
for a theory to describe this transition smoothly. Figure [T] shows the Born- 
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Fig. 1. Potential energy surface of C2 for various SIC calculations as indicated. 

Oppenheimer potential energy surface for the C2 dimer computed with a va- 
riety of approaches. Let us first start with the LDA approach which provides 
a qualitatively good approach with a fair reproduction of both bond length 
and dissociation energy but which unfortunately underestimates the ground 
state vibration frequency by about 25 %. All SIC corrected methods provide 
a much better reproduction of the bond length and the dissociation energy of 
the equilibrium state (less than 5 % of discrepancy with respect to the exper- 
imental data). They also improve the value of the vibration frequency : Slater 
and ADSIC yield an agreement within less than 5 %, and full SIC, GS(var) 
and GS(sym) within typically 10-15 %. Note also that GS(var) deliver re- 
sults which are almost identical to GS(sym) while being much less expensive 
numerically. This will also hold in the dynamical regime (see Sec. 14. 4p . 
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However one observes strange behaviors in ADSIC and Slater at intermediate 
distance with the appearance of a totally unphysical "bump" in the poten- 
tial energy surface. The effect is not present neither in full SIC nor in GS 
approximations which thus both provide a correct account of the potential 
energy surface. The defect observed in Slater and ADSIC has different origin. 
In Slater, it is probably to be attributed to a conflict between a tendency of 
the system to create " delocalized" orbitals, to ensure bonding, and a tendency 
towards "localized" orbitals, to ensure a better account of the SIC. The two set 
strategy proves here very valuable by resolving the conflict. The ADSIC prob- 
lem comes from the fact that asymptotically the ADSIC correction should 
take a form different from that at smaller distance because of the different 
number of involved electrons (4 in each separate C atom, 8 in the dimer). 
ADSIC requires compact systems and is generally not suited for describing 
fragmentation. 



4-3 Polarizabilities 



Polar izabilities are a sensitive test case for density functional approaches 
[36 | [591465] . We will thus discuss this issue here for three sufficiently differ- 
ent systems. As a first example, we consider the carbon molecules C2 and C4. 
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Fig. 2. Transverse and longitudinal polarizabilities of the C2 molecule (left) and 
the C4 chain (right), calculated in various SIC schemes. Horizontal lines emphasize 
the SIC benchmark values and ease the comparison with the other results. 

Figure [2] shows their polarizabilities for the various approximations. To put 
the subsequent results on C molecules into perspective, we recall the computed 
polarizations for the C atom: along z axis, Oz = 10.40 Oq^ for both SIC and 
GS, along x,y axes, a^^y = 11.52 qq^ for GS and 11.76 Gq^ for SIC. Experi- 
mental values for the molecular polarizabilities seem not be available. But one 
can compare with other computed values obtained with much different meth- 
ods [66l[67] They yield generally comparable values. In [66], the longitudinal 
polarizability for C2 is ay = 25 ag^ for the ab initio methods and 34 a^^ for 
LDA/GGA, while the transverse one is a± =25 Oq^ or 100 Oq^ respectively, 
the latter value being a strange exception. The results for C4 are a 11 = 92 or 
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94 ao^ and a± = 30 or 32 oq^. Our results are generally lower for a±. However, 
it is to be noted that our calculations differ in the employed functionals and 
pseudopotentials which both can have a sensitive influence on the results. In 
view of that, the comparison as a whole looks satisfying. 

The main aim of figure [2] is a comparison of methods within the same setup. 
The C2 dimer shows the larger variance of the results and is obviously more 
critical than the C4 chain. It is obvious that GS provides the best approxima- 
tions to full SIC and it is interesting that both versions, GS(sym) and GS(var), 
perform almost equally well. 
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Fig. 3. Polarizabilities, for various SIC schemes as indicated, of the ground state 
configuration H4 in the T-shaped configuration, displayed in the right panel. Hori- 
zontal lines emphasize the SIC benchmark values and ease the comparison with the 
other results. 

The H4 ground state configuration is a "T-shaped" molecule [68] as indicated 
in the right panel of figure [3l The H4 ground state configuration consists of 
two H2 dimers bound with a H2-II2 center of mass distance of 6.425 oq. This is 
a demanding configuration as it contains two well localized cloud of electrons 
at each H2 center loosely connected between the centers. Traditional SIC- 
Slater and KLI tends to delocalize the wave functions too much. The triaxial 
spatial configuration provides three different polarizabilities depending on the 
orientation of the external electric field relative to the molecule. The left panel 
of figure 13] compares the results for the three polarizabilities. Again we see that 
both variants of GS come very close to the benchmark (SIC), while LDA, as 
well as ADSIC, overestimate the polarizabilities, and traditional SIC-Slater is 
totally off. The overestimation is related to an exaggerated delocalization for 
LDA and ADSIC. The failure of the traditional Slater approximation indicates 
a too strong localization of the bonds. GS finds the right compromise. 
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Fig. 4. Polarizabilities, for various SIC schemes as indicated, of Nas, displayed in 
the right panel. Horizontal lines emphasize the SIC benchmark values and ease the 
comparison with the other results. 

As a final test case for polarizability, we consider the small sodium cluster 
Nas representative of simple metallic systems. Being a piece of a metal this 
cluster should have a delocalized electron cloud. But the rather soft binding 
of Nas degrades metalicity and drives towards weak localization. This makes 
Nas a particularly critical test case amongst metallic clusters [15] . The cluster 
is planar (see the right panel of figure Hj) which corresponds to a triaxial shape 
and leads to three rather different polarizabilities along the three major axes 
of the system. In order to display in a better readable way the various results, 
we have thus chosen to present values relative to the full SIC ones rather 
than absolute values. The left panel of figure H] shows the polarizabilities. 
We obtain much larger absolute values of polarizabilities than in the case of 
organic systems due to the metallic nature of bonding (delocalization and 
lower binding). Not surprisingly, all approaches perform rather well, better 
than in organic systems, as is to be expected for a simple metallic system. 
Within the lower error bands, we still see differences in the performance with 
a clear improvement provided by both GS versions. 



4-4 Time- dependent case 



As discussed in sections 13.2.31 and 13.2.41 Generalized Slater does not exactly 
fulfill conservation laws and has thus to be used with caution in actual time 
dependent processes. Still it is interesting to test in dynamical scenarios, espe- 
cially in the linear domain where it could advantageously replace more compli- 
cated approaches. The term linear domain refers to electronic oscillations with 
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small amplitudes. In the context of clusters and molecules, it largely refers to 
the analysis of the optical response which plays a key role in a broad vari- 
ety of dynamical scenarios, both in the linear and non linear domains [TT] . It 
thus represents a key issue in these systems. TDDFT, in particular in its real 
time formulation [69], is especially well suited to address such phenomena. For 
then, the point is simply to excite the system with a sufficiently small energy, 
whatever its value, so that dynamics sets in and allows a spectral analysis [69] . 
We shall thus mostly discuss such cases in the following. 



4.4-1 Na5 

As a first test case for dynamics, we consider Na5 which was found to be 
critical probe for studying conservation properties [15] because of its soft and 
easily polarizable electron cloud. We excite the electronic cloud by applying 
a boost in the x direction to each wave function. This simulates a very short 
laser pulse. We compare the case of very small excitation in the linear regime 
with that of a larger excitation. No absorbing boundary is used here, so that 
the total energy should be conserved in time. 

We first start with the low excitation case presented in Figure O The left panel 




30 60 90 120 150 180 30 60 90 120 150 180 
time (fs) time (fs) 

Fig. 5. (Color online) Time evolution of the total energy (left) and the dipole mo- 
ment in the spatial x direction (right) for Nas after a boost of the wave functions 
with momentum 0.001/ao in the x direction, for various SIC calculations as indi- 
cated. For the sake of clarity, some results have been down-shifted by a constant 
offset. 

compares energy conservation. In variational approaches, the total energy is 
conserved. We thus plot the deviation to this energy conservation, that is 
AE/E = [E{t)- E{0)]/E{0). As expected, LDA and TDSIC show up as 
straight lines because these methods are proven to conserve energy. TDGS 
will finally also develop an energy instability but it stays stable for much 
longer more than 7 times the standard Slater approximation which diverges 
after only 30 fs. Moreover, we see that the faster variant GS(var) performs 
as well as the version employing the involved symmetry condition. The right 
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panels of figure |5] show tlie time evolution of dipole moments. For the sake of 
clarity, the signals for GS(sym) and GS(var) have been shifted. In agreement 
with the time evolution of the total energy, in a TDGS calculation, the dipoles 
do not exhibit any significant evolution, while a standard Slater calculation 
produces large oscillations. 

Figure [6] shows the energy conservation for the same test case Na5 but for a 
much larger excitation energy 50 times higher, thus in the non-linear regime. 
We see again that both variants of TDGS provide a longer stability time 
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Fig. 6. (Color online) Total energy deviation of Nas for a boost momentum of 
0.05/ao in the x direction. 

than the standard Slater approximation. However, the stability time is much 
shorter than in the previous case of small excitation. The quality of Slater 
approximations is degrading with increasing excitation. The applicability has 
to be checked for each system and dynamical range anew. 



4-4-2 T-shape H^ 

We now turn to the case of H4 in the T configuration, as displayed in the 
top right panel of Figure [3 We excite this cluster the same way, that is with 
an instantaneous boost. Since GS(var) and GS(sym) turn out to give very 
close results, we only plot one of these results and denote them as GS without 
mentioning the criterion used for the calculation of the unitary transform. The 
findings from the lower two panels are the same as from the previous results. 
The standard Slater approximation runs rather quickly into violation of energy 
conservation, while GS stays stable for much longer. The instability in energy 
shows up at later times in the dipole signal. The upper left panel explores 
energy conservation for GS at various initial excitation energies. It is apparent 
that time span at which GS propagates in stable manner depends crucially on 
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Fig. 7. (Color online) T-shaped H4 (ionic configuration in the top right panel) 
excited by boosts bx in the x direction as indicated. Left column: total energy 
deviation as a function of time; bottom right panel: time evolution of the x dipole 
moment. 

the excitation energy of the process. At the present stage of development it 
can be safely used only in the linear domain. 



5 Conclusion 



We have presented in this paper a generalized formulation of SIC-OEP in the 
time domain. It relies on the introduction of a double set of electron orbitals. 
The double-set strategy plays a key role for the time propagation of full SIC 
but also provides in the OEP context a valuable tool for deriving approximate 
versions of the full SIC-OEP, in stationary as well as time dependent processes. 
It allows in particular to introduce the Generalized Slater (GS) approximation 
which preserves the simplicity of the standard Slater approximation and yet 
extends its range of validity. While a first set of " in general localized" orbitals 
serves to fulfill ortho-normality and to construct the approximate Hamiltonian, 
a second set of "physical" wavefunctions allows a simple time propagation. We 
have shown that the GS improves over the defects of the simple Slater approx- 
imation in particular what concerns conservation laws but still does not allow 
to fulfill them exactly. We have thus performed various tests on a series of 
clusters and molecules in order to explore in a practical way the capabilities 
of GS. We have seen that it performs quite well for static properties and al- 
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lows to deal with dynamics in the linear domain. It then provides a simple and 
transparent alternative to the full TD-OEP or TD-SIC which require substan- 
tial numerical effort. The numerical performance of the method still require 
to be optimized. We have shown that the expensive symmetry condition may 
be advantageously replaced by the conceptually and computationally simpler 
"localization" ansatz which opens up possibly new, simplified, approaches to 
the problem. This strategy nevertheless still requires to be explored in more 
detail. Work along that line is in progress. 
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